% function [FEAS_FLAG, bu_a, info] = minPICU(H, bu_min, w, cb, opts)
% minPICU
%
% U-user / 1-tone / scalar Gaussian interference-channel prototype.
%
% This routine searches over full subset IC atomizations and receiver
% decodings.
%
% -------------------------------------------------------------------------
% INPUTS
% -------------------------------------------------------------------------
%
% H
%   U x U scalar channel matrix.
%
%   H(r,u) is the scalar channel gain from transmitter/user u to receiver r.
%
% bu_min
%   U x 1 or 1 x U vector of requested user rates, in bits per channel use.
%
%   bu_min(u) is the requested total rate for user u, summed over all atoms
%   belonging to user u.
%
% w
%   U x 1 or 1 x U positive weight vector for user energies.
%
%   The optimized weighted energy is
%
%       E_total = sum_u w(u) * E_users(u).
%
% cb
%   Baseband convention:
%
%       cb = 1  complex baseband, rate = log2(1 + SINR)
%       cb = 2  real baseband,    rate = 0.5*log2(1 + SINR)
%
% opts
%   Optional structure controlling tolerances, iteration counts, verbosity,
%   order-search limits, time-sharing, and ellipsoid settings.
%
% -------------------------------------------------------------------------
% OUTPUTS
% -------------------------------------------------------------------------
%
% FEAS_FLAG
%   Feasibility/status flag:
%
%       0  target was not met by the harvested vertex/time-sharing pool
%       1  target was met by a single vertex
%       2  target was met by time-sharing among harvested vertices
%
% bu_a
%   1 x U achieved user-rate vector after final selection/time-sharing.
%
%   bu_a(u) is the achieved total rate for user u.
%
% info
%   Main output structure. Important fields are:
%
%   info.mode
%       Text label describing the final solution type, for example
%       'single-vertex', 'time-share-lp', or 'target-miss'.
%
%   info.target_requested
%       U x 1 requested user-rate vector.
%
%   info.bu_a
%       U x 1 achieved user-rate vector.
%
%   info.E
%       Final weighted total energy:
%
%           info.E = sum_u w(u)*info.E_users(u).
%
%   info.E_users
%       U x 1 vector of total transmit energy per user.
%
%       info.E_users(u) is the sum of the powers of all atoms belonging to
%       user u.
%
%   info.powers
%       K x 1 vector of atom/subuser powers for the final averaged solution.
%
%       K = U^2 in this pairwise IC atomization.
%
%   info.atom_energies
%       Alias for info.powers. This is often the clearer name when passing
%       the result to receiver-design code.
%
%   info.atom_rates
%       K x 1 vector of reliable atom/subuser rates for the final averaged
%       solution.
%
%       For an atom decoded at more than one receiver, the atom rate is the
%       minimum over those receivers.
%
%   info.subuser_rates
%       Alias for info.atom_rates.
%
%   info.atom_rates_support
%       K x 1 vector of support-function atom rates. In the current exact
%       Imin version, this is forced equal to info.atom_rates.
%
%   info.atom_table
%       K x 1 structure array describing the atoms. For atom k:
%
%           info.atom_table(k).index       atom index k
%           info.atom_table(k).label       readable label, e.g. 'u2_to_3'
%           info.atom_table(k).user        owning transmitter/user
%           info.atom_table(k).decoded_at  receivers that decode this atom
%           info.atom_table(k).type        'priv' or 'pair'
%           info.atom_table(k).home_rx     intended/home receiver
%           info.atom_table(k).mate_rx     paired receiver
%
%   info.schedule
%       Structure array describing the actual receiver-design schedules.
%
%       This is the most useful output for GDFE construction. If the final
%       point is a single vertex, numel(info.schedule)=1. If the final point
%       is a time-sharing solution, each entry is one active time-sharing
%       vertex.
%
%       For schedule entry s:
%
%           info.schedule(s).lambda
%               Time-sharing fraction for this vertex.
%
%           info.schedule(s).atom_energies
%               K x 1 atom powers used in this vertex.
%
%           info.schedule(s).atom_rates
%               K x 1 atom rates used in this vertex.
%
%           info.schedule(s).bu
%               U x 1 user-rate vector produced by this vertex.
%
%           info.schedule(s).E
%               Weighted total energy for this vertex.
%
%           info.schedule(s).E_users
%               U x 1 per-user energies for this vertex.
%
%           info.schedule(s).rx_orders
%               U x 1 cell array. rx_orders{r} is the ordered list of atom
%               indices decoded at receiver r.
%
%           info.schedule(s).orders
%               Same order information in the internal orders structure.
%
%           info.schedule(s).rx_rate
%               K x U matrix. rx_rate(k,r) is the rate assigned to atom k
%               when decoded at receiver r. Entries are zero when atom k is
%               not decoded at receiver r.
%
%           info.schedule(s).active_rx
%               K x 1 cell array. active_rx{k} lists the receiver(s) that
%               attain the Imin minimum for atom k.
%
%   info.theta_final
%       Final ellipsoid multiplier vector.
%
%   info.theta_history
%       U x N matrix of multiplier iterates.
%
%   info.ellip_trace
%       Iteration-by-iteration diagnostic trace.
%
%   info.vertex_pool_*
%       Raw and lower-envelope harvested candidate vertices.
%
%   info.face_pool_*
%       Extra tied-theta face vertices, when face mode is active.
%
%   info.lp_lambda, info.lp_active_idx
%       LP time-sharing fractions and active vertex indices, if time-sharing
%       was used.
%
% -------------------------------------------------------------------------
% ATOMIZATION
% -------------------------------------------------------------------------
%
% The program uses the full IC subset atomization:
%
%   For each user u, there is one atom for every receiver subset S
%   such that u is in S.
%
% Number of atoms per user:
%
%       2^(U-1)
%
% Total number of atoms:
%
%       K = U * 2^(U-1)
%
% For U = 3:
%
%       K = 3 * 4 = 12 atoms
%
% At receiver r, the decoded atoms are exactly those whose decoded_at
% set contains r. All other active atoms are treated as noise.
%
% -------------------------------------------------------------------------
function [FEAS_FLAG, bu_a, info] = minPICU(H, bu_min, w, cb, opts)
    if nargin < 5
        opts = struct();
    end
    opts = mp_set_defaults(opts);

    if isfield(opts, 'randomSeed')
        rng(opts.randomSeed);
    end

    [Ur, Uc] = size(H);
    assert(Ur == Uc, 'H must be square.');
    U = Ur;

    bu_min = bu_min(:);
    w = w(:);

    assert(numel(bu_min) == U, 'bu_min must have length U.');
    assert(numel(w) == U && all(w > 0), 'w must have length U and be positive.');
    assert(cb == 1 || cb == 2, 'cb must be 1 or 2.');

    G = abs(H).^2;
    atoms = mp_build_atom_table(U);
    K = numel(atoms);

    if opts.verbose
        fprintf('minPICU:\n');
        fprintf('  U = %d, K = %d atoms\n', U, K);
        fprintf('  target = [%s]\n', num2str(bu_min.', ' %.6g'));
        fprintf('  weights = [%s]\n', num2str(w.', ' %.6g'));
        fprintf('  cb = %d\n', cb);
        fprintf('  exact Imin reward used everywhere = 1\n');
    end

    theta_max = mp_compute_theta_max_bounds(G, bu_min, w, cb, opts, atoms);
    theta_max = max(theta_max(:), opts.thetaFloor);
    theta0 = 0.5 * theta_max(:);

    xc = theta0(:);
    A = diag(theta_max(:).^2);

    if opts.verbose
        fprintf('  theta_max = [%s]\n', num2str(theta_max.', ' %.6g'));
        fprintf('  theta0    = [%s]\n', num2str(theta0.', ' %.6g'));
    end

    theta_history = zeros(U, opts.maxEllipIter);
    ellip_trace = repmat(struct( ...
        'iter', [], ...
        'theta', [], ...
        'nOrders', [], ...
        'best_bu', [], ...
        'best_bu_support', [], ...
        'best_E', [], ...
        'shortfall', [], ...
        'pool_idx', []), opts.maxEllipIter, 1);

    vertex_pool = mp_init_vertex_pool(K, U);
    prev_best = [];

    %----------------------------------------------------------------------
    % Main ellipsoid loop (exact-rate driven)
    %----------------------------------------------------------------------

    for it = 1:opts.maxEllipIter
        theta_history(:, it) = xc;

        orderSet = mp_generate_order_candidates(xc, atoms, U, opts, prev_best);

        bestIter = [];
        idx_kept = [];

        for kk = 1:numel(orderSet)
            cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orderSet{kk}, xc, prev_best);

            [vertex_pool, idxAdded] = mp_pool_add_vertex(vertex_pool, cand, opts.poolUniqueTol);
            idx_kept = [idx_kept; idxAdded]; %#ok<AGROW>

            if isempty(bestIter) || cand.lag_obj < bestIter.lag_obj
                bestIter = cand;
            end
        end

                prev_best = bestIter;
        shortfall = bu_min(:) - bestIter.bu(:);

        %--------------------------------------------------------------
        % Adaptive theta-box expansion:
        % if a user still has meaningful shortfall and its multiplier is
        % pressed near the current theta_max wall, enlarge that wall.
        %--------------------------------------------------------------
        if opts.thetaAdaptiveExpand
            hit = (xc >= opts.thetaNearMaxFrac * theta_max) & ...
                  (shortfall > opts.thetaExpandShortTol);

            if any(hit)
                old_theta_max = theta_max;
                theta_max(hit) = min(opts.thetaMaxCap, ...
                                     opts.thetaExpandFactor * theta_max(hit));

                % Preserve the ellipsoid scaling when the box expands.
                for uu = 1:numel(theta_max)
                    if theta_max(uu) > old_theta_max(uu)
                        scaleu = theta_max(uu) / max(old_theta_max(uu), opts.thetaFloor);
                        A(uu,uu) = max(A(uu,uu), theta_max(uu)^2);
                        xc(uu) = min(xc(uu), theta_max(uu));
                    end
                end

                if opts.verbose
                    fprintf('    expanded theta_max -> [%s]\n', ...
                        num2str(theta_max.', ' %.6g'));
                end
            end
        end

        ellip_trace(it).iter = it;
        ellip_trace(it).theta = xc(:);
        ellip_trace(it).nOrders = numel(orderSet);
        ellip_trace(it).best_bu = bestIter.bu(:);
        ellip_trace(it).best_bu_support = bestIter.bu_support(:);
        ellip_trace(it).best_E = bestIter.E;
        ellip_trace(it).shortfall = shortfall(:);
        ellip_trace(it).pool_idx = idx_kept(:).';

        if opts.verbose
            fprintf(['  iter %2d: theta=[%s], orders=%d, E=%-.6g, ' ...
                     'bu=[%s], bu_sup=[%s], short=[%s]\n'], ...
                it, num2str(xc.', ' %.5g'), numel(orderSet), bestIter.E, ...
                num2str(bestIter.bu.', ' %.5g'), ...
                num2str(bestIter.bu_support.', ' %.5g'), ...
                num2str(shortfall.', ' %.3e'));
        end

        % Stop on exact-rate hit.
        if norm(shortfall, inf) <= opts.hitTol
            break;
        end

        g = shortfall(:);
        if norm(g) <= 1e-14
            break;
        end

        denom = sqrt(g.' * A * g);
        if denom <= 1e-14
            break;
        end

        n = U;
        Ag = A * g;
        xc = xc + (1/(n+1)) * (Ag / denom);
        A = (n^2/(n^2-1)) * (A - (2/(n+1)) * ((Ag * Ag.') / max(g.' * A * g, 1e-30)));

        xc = min(max(xc, 0), theta_max);
        A = mp_symmetrize_pd(A, theta_max);
    end

    theta_history = theta_history(:, 1:it);
    ellip_trace = ellip_trace(1:it);
    %----------------------------------------------------------------------
    % Extra exploratory harvesting:
    % The LP can only succeed if the harvested vertex pool is rich enough.
    % For U>=3, add several perturbed-theta sweeps to expose additional
    % corner/face vertices before final convexification.
    %----------------------------------------------------------------------
    if U >= 3
        thetaSeeds = {};
        thetaSeeds{end+1} = xc(:);
        thetaSeeds{end+1} = 0.75 * xc(:);
        thetaSeeds{end+1} = 1.25 * xc(:);

        for uu = 1:U
            e = zeros(U,1);
            e(uu) = 1;
            thetaSeeds{end+1} = xc(:) + 0.25 * max(1, norm(xc,inf)) * e;
        end

        for uu = 1:U
            e = zeros(U,1);
            e(uu) = 1;
            thetaSeeds{end+1} = max(0, xc(:) - 0.25 * max(1, norm(xc,inf)) * e);
        end

        for ss = 1:numel(thetaSeeds)
            th = min(max(thetaSeeds{ss}, 0), theta_max);
            orderSet_extra = mp_generate_order_candidates(th, atoms, U, opts, prev_best);

            for kk = 1:numel(orderSet_extra)
                cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orderSet_extra{kk}, th, prev_best);
                [vertex_pool, ~] = mp_pool_add_vertex(vertex_pool, cand, opts.poolUniqueTol);
            end
        end
    end
    %----------------------------------------------------------------------
    % Final selection: if there are tied theta clusters, harvest tied face;
    % otherwise use global select over harvested CIC vertex pool.
    %----------------------------------------------------------------------

    used_face_mode = false;
    face_pool = mp_init_vertex_pool(K, U);

    if opts.enableFaceMode && mp_has_tied_theta_cluster(xc, opts.tieTol)
        face_pool = mp_harvest_tied_theta_face(G, bu_min, w, cb, opts, atoms, xc, prev_best);
        face_pool = mp_pool_lower_envelope(face_pool, opts);
        [sel_face, ok_face] = mp_global_select(face_pool, bu_min, opts.hitTol, opts.degenerateTol, opts);
        if ok_face > 0
            sel = sel_face;
            FEAS_FLAG = ok_face;
            if FEAS_FLAG == 1
                sel.mode = 'single-vertex-face';
            else
                sel.mode = 'time-share-face';
            end
            used_face_mode = true;
        end
    end

    raw_vertex_pool = vertex_pool;

    if ~used_face_mode
      [sel, FEAS_FLAG] = mp_global_select(raw_vertex_pool, bu_min, opts.hitTol, opts.degenerateTol, opts);
    end

vertex_pool = raw_vertex_pool;   % keep raw pool for diagnostics and schedule safety

    bu_a = sel.bu(:).';

    info = struct();
    info.mode = sel.mode;
    info.target_requested = bu_min(:);
    info.target_supported = mp_tern(FEAS_FLAG > 0, bu_min(:), []);
    info.bu_a = bu_a(:);

    info.E = sel.E;
    info.E_users = sel.E_users(:);
    info.powers = sel.powers(:);
    info.atom_rates = sel.atom_rates(:);
    info.atom_rates_support = sel.atom_rates_support(:);
    info.atom_table = atoms;

    info.theta_max = theta_max(:);
    info.theta0 = theta0(:);
    info.theta_final = xc(:);
    info.theta_history = theta_history;
    info.ellip_trace = ellip_trace;

    info.vertex_pool_bu = vertex_pool.bu;
    info.vertex_pool_bu_support = vertex_pool.bu_support;
    info.vertex_pool_E = vertex_pool.E(:);
    info.vertex_pool_meta = {vertex_pool.meta{:}};
    pruned_pool = mp_pool_lower_envelope(vertex_pool, opts);
    info.vertex_pool_lower_bu = pruned_pool.bu;
    info.vertex_pool_lower_E = pruned_pool.E(:);
    info.vertex_pool_lower_meta = {pruned_pool.meta{:}};

    info.face_pool_bu = face_pool.bu;
    info.face_pool_bu_support = face_pool.bu_support;
    info.face_pool_E = face_pool.E(:);
    info.face_pool_meta = {face_pool.meta{:}};

    info.best_vertex = sel.best_vertex;
    info.ts_exact = sel.ts_exact;
    info.ts_feasible = sel.ts_feasible;
    info.lp_lambda = sel.lp_lambda;
    info.lp_active_idx = sel.lp_active_idx;
    info.lp_exitflag = sel.lp_exitflag;
    info.lp_message  = sel.lp_message;
    info.G = G;
    info.U = U;
    info.K = K;
    info.used_face_mode = used_face_mode;

    % ---------------------------------------------------------------------
    % Receiver-design export fields.
    % These aliases/schedules do not affect optimization; they only expose
    % final outputs for downstream GDFE/SIC construction.
    % ---------------------------------------------------------------------
    info.atom_energies = info.powers;        % K x 1 final averaged atom powers
    info.subuser_rates = info.atom_rates;    % K x 1 final averaged atom rates

    info.schedule = mp_build_schedule_from_selection(sel, vertex_pool, face_pool, used_face_mode);

    if opts.verbose
        fprintf('\nFinal result:\n');
        fprintf('  FEAS_FLAG = %d\n', FEAS_FLAG);
        fprintf('  bu_a      = [%s]\n', num2str(bu_a, ' %.6g'));
        fprintf('  mode      = %s\n', info.mode);
        fprintf('  E_total   = %.6g\n', info.E);
        fprintf('  E_users   = [%s]\n', num2str(info.E_users.', ' %.6g'));
        fprintf('  theta     = [%s]\n', num2str(info.theta_final.', ' %.6g'));
        if info.ts_exact.used
            fprintf('  ts_exact.used = 1\n');
            fprintf('  ts_exact.frac = [%s]\n', num2str(info.ts_exact.frac(:).', ' %.6g'));
            fprintf('  ts_exact.idx  = [%s]\n', num2str(info.ts_exact.idx(:).', ' %d'));
        elseif info.ts_feasible.used
            fprintf('  ts_feasible.used = 1\n');
            fprintf('  ts_feasible.frac = [%s]\n', num2str(info.ts_feasible.frac(:).', ' %.6g'));
            fprintf('  ts_feasible.idx  = [%s]\n', num2str(info.ts_feasible.idx(:).', ' %d'));
        elseif ~isempty(info.lp_lambda)
            fprintf('  lp time-sharing used\n');
            fprintf('  lp active idx = [%s]\n', num2str(info.lp_active_idx(:).', ' %d'));
            fprintf('  lp lambda     = [%s]\n', num2str(info.lp_lambda(:).', ' %.6g'));
        else
            fprintf('  no time-sharing used\n');
        end
    end
end

%==========================================================================
% DEFAULTS / ATOMS / BASIC TABLES
%==========================================================================

function opts = mp_set_defaults(opts)
    if ~isfield(opts, 'verbose'), opts.verbose = 1; end
    if ~isfield(opts, 'maxEllipIter'), opts.maxEllipIter = 80; end
    if ~isfield(opts, 'hitTol'), opts.hitTol = 1e-6; end
    if ~isfield(opts, 'tieTol'), opts.tieTol = 5e-3; end
    if ~isfield(opts, 'poolUniqueTol'), opts.poolUniqueTol = 1e-8; end
    if ~isfield(opts, 'degenerateTol'), opts.degenerateTol = 1e-7; end
    if ~isfield(opts, 'nInnerStarts'), opts.nInnerStarts = 8; end
    if ~isfield(opts, 'innerMaxIter'), opts.innerMaxIter = 3000; end
    if ~isfield(opts, 'innerMaxFunEvals'), opts.innerMaxFunEvals = 12000; end
    if ~isfield(opts, 'innerTolX'), opts.innerTolX = 1e-8; end
    if ~isfield(opts, 'innerTolFun'), opts.innerTolFun = 1e-8; end
    if ~isfield(opts, 'bumpBits'), opts.bumpBits = 1; end
    if ~isfield(opts, 'thetaInflation'), opts.thetaInflation = 10.0; end
    if ~isfield(opts, 'thetaFloor'), opts.thetaFloor = 1e-8; end
    if ~isfield(opts, 'thetaAdaptiveExpand'), opts.thetaAdaptiveExpand = true; end
    if ~isfield(opts, 'thetaExpandFactor'), opts.thetaExpandFactor = 2.0; end
    if ~isfield(opts, 'thetaNearMaxFrac'), opts.thetaNearMaxFrac = 0.98; end
    if ~isfield(opts, 'thetaExpandShortTol'), opts.thetaExpandShortTol = 1e-3; end
    if ~isfield(opts, 'thetaMaxCap'), opts.thetaMaxCap = 1e6; end    
    if ~isfield(opts, 'capAtTarget'), opts.capAtTarget = true; end
    if ~isfield(opts, 'enableFaceMode'), opts.enableFaceMode = true; end
    if ~isfield(opts, 'faceUniqueTol'), opts.faceUniqueTol = 1e-7; end
    if ~isfield(opts, 'useAvgSupportForUge3'), opts.useAvgSupportForUge3 = false; end
    if ~isfield(opts, 'maxTiePermsPerRx'), opts.maxTiePermsPerRx = 24; end
    if ~isfield(opts, 'maxGlobalOrderCombos'), opts.maxGlobalOrderCombos = 200; end
    if ~isfield(opts, 'nInnerStarts'), opts.nInnerStarts = 20; end
    if ~isfield(opts, 'orderScoreTol'), opts.orderScoreTol = 1e-6; end
    if ~isfield(opts, 'nRandomOrderSweeps'), opts.nRandomOrderSweeps = 12; end
    if ~isfield(opts, 'orderScoreEps'), opts.orderScoreEps = 1e-8; end
    if ~isfield(opts, 'orderPrevMix'), opts.orderPrevMix = 0.5; end
    if ~isfield(opts, 'lpTol'), opts.lpTol = 1e-9; end
    if ~isfield(opts, 'poolDomTol'), opts.poolDomTol = 1e-9; end
    if ~isfield(opts, 'enableHullPrune'), opts.enableHullPrune = true; end
    if ~isfield(opts, 'hullPruneEvery'), opts.hullPruneEvery = 25; end
    if ~isfield(opts, 'compressActiveSet'), opts.compressActiveSet = true; end
end

function atoms = mp_build_atom_table(U)
% Full IC subset atomization.
%
% Atom for user u is decoded by every receiver subset S such that u in S.
%
% Number of atoms per user:
%     2^(U-1)
%
% Total atoms:
%     K = U * 2^(U-1)
%
% Number decoded at any receiver r:
%     own-user atoms:          2^(U-1)
%     other-user atoms: (U-1)*2^(U-2)
%     total:            (U+1)*2^(U-2)
%
% For U=3:
%     K = 3*4 = 12 total atoms
%     each receiver decodes 8 atoms max

    K = U * 2^(U-1);

    atoms = repmat(struct( ...
        'index', [], ...
        'label', '', ...
        'user', [], ...
        'decoded_at', [], ...
        'home_rx', [], ...
        'mate_rx', [], ...
        'type', ''), K, 1);

    k = 0;

    for u = 1:U

        others = setdiff(1:U, u);

        % Enumerate all subsets containing u.
        % mask selects which "other" receivers also decode user u's atom.
        for mask = 0:(2^(U-1)-1)

            extra = [];
            for j = 1:numel(others)
                if bitget(mask, j)
                    extra(end+1) = others(j); %#ok<AGROW>
                end
            end

            decset = sort([u extra]);

            k = k + 1;
            atoms(k).index = k;
            atoms(k).user = u;
            atoms(k).home_rx = u;
            atoms(k).decoded_at = decset;
            atoms(k).mate_rx = setdiff(decset, u);

            if numel(decset) == 1
                atoms(k).type = 'priv';
                atoms(k).label = sprintf('u%d_priv', u);
            elseif numel(decset) == 2
                atoms(k).type = 'pair';
                atoms(k).label = sprintf('u%d_to_%d', u, atoms(k).mate_rx);
            else
                atoms(k).type = 'subset';
                atoms(k).label = sprintf('u%d_to_%s', u, sprintf('%d', atoms(k).mate_rx));
            end
        end
    end
end

function idx = mp_user_atom_indices(atoms, u)
    idx = find([atoms.user] == u);
end

function idx = mp_decodable_atoms_at_receiver(atoms, r)
    K = numel(atoms);
    tf = false(K,1);
    for k = 1:K
        tf(k) = any(atoms(k).decoded_at == r);
    end
    idx = find(tf).';
end

function pool = mp_init_vertex_pool(K, U)
    pool = struct();
    pool.bu = zeros(0, U);
    pool.bu_support = zeros(0, U);
    pool.E = zeros(0, 1);
    pool.E_users = zeros(0, U);
    pool.powers = zeros(0, K);
    pool.atom_rates = zeros(0, K);
    pool.atom_rates_support = zeros(0, K);
    pool.meta = {};
    pool.n_add = 0;
end


function [pool, idx] = mp_pool_add_vertex(pool, ev, tolU)
    idx = [];

    if isempty(pool.bu)
        pool = mp_pool_append_vertex(pool, ev);
        idx = 1;
        return;
    end

    % Keep only the lowest-energy realization for essentially identical b.
    d = sqrt(sum((pool.bu - ev.bu(:).').^2, 2));
    same = find(d <= tolU, 1, 'first');

    if ~isempty(same)
        if ev.E + 1e-14 < pool.E(same)
            pool = mp_pool_replace_vertex(pool, same, ev);
        end
        idx = same;
        return;
    end

    % Do not prune the raw harvested pool online. The previous revision
    % became too aggressive for U>2 and removed the very diversity needed
    % for the final time-sharing LP.
    pool = mp_pool_append_vertex(pool, ev);
    idx = size(pool.bu,1);
end

function pool = mp_pool_append_vertex(pool, ev)
    pool.bu(end+1,:) = ev.bu(:).';
    pool.bu_support(end+1,:) = ev.bu_support(:).';
    pool.E(end+1,1) = ev.E;
    pool.E_users(end+1,:) = ev.E_users(:).';
    pool.powers(end+1,:) = ev.powers(:).';
    pool.atom_rates(end+1,:) = ev.atom_rates(:).';
    pool.atom_rates_support(end+1,:) = ev.atom_rates_support(:).';
    pool.meta{end+1,1} = mp_make_meta(ev);
end

function pool = mp_pool_replace_vertex(pool, i, ev)
    pool.bu(i,:) = ev.bu(:).';
    pool.bu_support(i,:) = ev.bu_support(:).';
    pool.E(i) = ev.E;
    pool.E_users(i,:) = ev.E_users(:).';
    pool.powers(i,:) = ev.powers(:).';
    pool.atom_rates(i,:) = ev.atom_rates(:).';
    pool.atom_rates_support(i,:) = ev.atom_rates_support(:).';
    pool.meta{i} = mp_make_meta(ev);
end

function tf = mp_point_dominates(b1, E1, b2, E2, tol)
    better_or_equal_rates = all(b1(:) >= b2(:) - tol);
    no_worse_energy = (E1 <= E2 + tol);
    strictly_better_somewhere = any(b1(:) > b2(:) + tol) || (E1 < E2 - tol);
    tf = better_or_equal_rates && no_worse_energy && strictly_better_somewhere;
end

function pool = mp_pool_filter(pool, keep)
    keep = keep(:);
    pool.bu = pool.bu(keep,:);
    pool.bu_support = pool.bu_support(keep,:);
    pool.E = pool.E(keep,:);
    pool.E_users = pool.E_users(keep,:);
    pool.powers = pool.powers(keep,:);
    pool.atom_rates = pool.atom_rates(keep,:);
    pool.atom_rates_support = pool.atom_rates_support(keep,:);
    pool.meta = pool.meta(keep,:);
end

function pool = mp_pool_lower_envelope(pool, opts)
    if isempty(pool.bu) || size(pool.bu,1) <= 2
        return;
    end
    if ~isfield(opts,'enableHullPrune') || ~opts.enableHullPrune
        return;
    end
    if ~license('test','optimization_toolbox')
        return;
    end

    nV = size(pool.bu,1);
    keep = true(nV,1);
    tol = max(getfield(opts,'lpTol'), 1e-9);

    % A point i is discardable if the remaining points can achieve
    % b >= b_i with weighted energy <= E_i.
    for i = 1:nV
        idx = find(keep);
        idx(idx == i) = [];
        if isempty(idx)
            continue;
        end

        B = pool.bu(idx,:);
        Evec = pool.E(idx);
        target_i = pool.bu(i,:).';
        Ei = pool.E(i);

        f = Evec(:);
        Aineq = -B.';
        bineq = -target_i;
        Aeq = ones(1, numel(idx));
        beq = 1;
        lb = zeros(numel(idx),1);
        ub = ones(numel(idx),1);

        try
            lpOpts = optimoptions('linprog', ...
                'Display','none', ...
                'OptimalityTolerance', tol, ...
                'ConstraintTolerance', tol);
            [lam, val, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, lpOpts);
            if exitflag > 0
                bi = B.' * lam;
                if all(bi >= target_i - 10*tol) && val <= Ei + 10*tol
                    keep(i) = false;
                end
            end
        catch
        end
    end

    pool = mp_pool_filter(pool, keep);
end

function val = mp_getfield_local(s, fld, defaultVal)
    if isstruct(s) && isfield(s, fld)
        val = s.(fld);
    else
        val = defaultVal;
    end
end

function meta = mp_make_meta(ev)
    meta = struct();
    meta.orders = ev.orders;
    meta.powers = ev.powers(:);
    meta.atom_rates = ev.atom_rates(:);
    meta.atom_rates_support = ev.atom_rates_support(:);
    meta.bu = ev.bu(:);
    meta.bu_support = ev.bu_support(:);
    meta.E = ev.E;
    meta.E_users = ev.E_users(:);
    meta.active_rx = ev.active_rx;
    meta.rx_rate = ev.rx_rate;
end

%==========================================================================
% THETA BOUNDS / STARTS / INNER SOLVES
%==========================================================================

function theta_max = mp_compute_theta_max_bounds(G, bu_min, w, cb, opts, atoms)
    U = numel(bu_min);
    theta_max = zeros(U,1);

    theta_probe = ones(U,1);
    probeOrders = mp_generate_order_candidates(theta_probe, atoms, U, opts, []);

    for u = 1:U
        bt = bu_min;
        bt(u) = bt(u) + opts.bumpBits;

        bestE = inf;
        for kk = 1:numel(probeOrders)
            cand = mp_inner_solve_targeted(G, bt, bu_min, w, cb, opts, atoms, probeOrders{kk});
            if cand.E < bestE
                bestE = cand.E;
            end
        end

        theta_max(u) = opts.thetaInflation * max(bestE, opts.thetaFloor);
    end
end

function starts = mp_build_inner_starts(G, bu_min, cb, opts, atoms, prev)
    U = numel(bu_min);
    K = numel(atoms);

    starts = zeros(K, opts.nInnerStarts);
    Pseed = zeros(K,1);

    for u = 1:U
        au = 2^(cb * bu_min(u)) - 1;
        p_iso = max(au / max(G(u,u), eps), 1e-10);
        idxu = mp_user_atom_indices(atoms, u);
        Pseed(idxu) = p_iso / numel(idxu);
    end

    starts(:,1) = log(max(Pseed, 1e-12));

    kk = 2;
    if ~isempty(prev) && kk <= opts.nInnerStarts
        starts(:,kk) = log(max(prev.powers(:), 1e-12));
        kk = kk + 1;
    end

    while kk <= opts.nInnerStarts
        scale = 10.^(0.8 * randn(K,1));
        starts(:,kk) = log(max(Pseed .* scale, 1e-12));
        kk = kk + 1;
    end
end

function cand = mp_inner_solve_targeted(G, bt, bu_min, w, cb, opts, atoms, orders)
    starts = mp_build_inner_starts(G, bu_min, cb, opts, atoms, []);
    best = [];

    for ss = 1:size(starts,2)
        z0 = starts(:,ss);
        f = @(zz) mp_targeted_objective(zz, G, bt, w, cb, opts, atoms, orders);

        z = fminsearch(f, z0, optimset( ...
            'Display', 'off', ...
            'TolX', opts.innerTolX, ...
            'TolFun', opts.innerTolFun, ...
            'MaxIter', opts.innerMaxIter, ...
            'MaxFunEvals', opts.innerMaxFunEvals));

        P = exp(z(:));
        ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);
        miss = max(bt(:) - ev.bu(:), 0);
        score = ev.E + 1e6 * sum(miss.^2);

        if isempty(best) || score < best.score
            best = ev;
            best.score = score;
        end
    end

    cand = best;
end

function val = mp_targeted_objective(z, G, bt, w, cb, opts, atoms, orders)
    P = exp(z(:));
    ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);
    miss = max(bt(:) - ev.bu(:), 0);
    val = ev.E + 1e6 * sum(miss.^2) + 1e-12 * sum(P.^2);
end

function cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orders, theta, prev)
    starts = mp_build_inner_starts(G, bu_min, cb, opts, atoms, prev);
    best = [];

    for ss = 1:size(starts,2)
        z0 = starts(:,ss);
        f = @(zz) mp_lagrangian_objective(zz, theta, G, bu_min, w, cb, opts, atoms, orders);

        z = fminsearch(f, z0, optimset( ...
            'Display', 'off', ...
            'TolX', opts.innerTolX, ...
            'TolFun', opts.innerTolFun, ...
            'MaxIter', opts.innerMaxIter, ...
            'MaxFunEvals', opts.innerMaxFunEvals));

        P = exp(z(:));
        ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);

        if opts.capAtTarget
            reward = min(ev.bu, bu_min);
        else
            reward = ev.bu;
        end

        ev.lag_obj = ev.E - theta(:).' * reward(:);

        if isempty(best) || ev.lag_obj < best.lag_obj
            best = ev;
        end
    end

    cand = best;
end

function val = mp_lagrangian_objective(z, theta, G, bu_min, w, cb, opts, atoms, orders)
    P = exp(z(:));
    ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders);

    if opts.capAtTarget
        reward = min(ev.bu, bu_min);
    else
        reward = ev.bu;
    end

    val = ev.E - theta(:).' * reward(:) + 1e-12 * sum(P.^2);
end

%==========================================================================
% RATE EVALUATION
%==========================================================================

function ev = mp_evaluate_fixed_orders(P, G, w, cb, opts, atoms, orders)
    K = numel(atoms);
    U = size(G,1);
    Cfun = @(x) (1/cb) * log2(1 + max(x,0));

    rx_rate = zeros(K, U);

    for r = 1:U
        ord = orders.rx_orders{r};
        done = false(1, K);
        decmask = false(1, K);
        decmask(ord) = true;

        for m = 1:numel(ord)
            k = ord(m);
            futureDec = decmask & ~done;
            futureDec(k) = false;

            denom = 1;
            for j = 1:K
                uj = atoms(j).user;
                if (~decmask(j)) || futureDec(j)
                    denom = denom + G(r, uj) * P(j);
                end
            end

            num = G(r, atoms(k).user) * P(k);
            rx_rate(k, r) = Cfun(num / denom);
            done(k) = true;
        end
    end

    atom_rates = zeros(K,1);
    atom_rates_support = zeros(K,1);
    bu = zeros(U,1);
    bu_support = zeros(U,1);
    active_rx = cell(K,1);

    for k = 1:K
        decset = atoms(k).decoded_at;
        rr = rx_rate(k, decset).';

        % Exact reliable rate
        rmin = min(rr);
        atom_rates(k) = rmin;

        % Support reward is forced equal to exact reliable rate
        atom_rates_support(k) = rmin;

        % Record active minimizing receivers
        act = decset(abs(rr - rmin) <= 1e-10);
        if isempty(act)
            act = decset(1);
        end
        active_rx{k} = act(:).';

        u = atoms(k).user;
        bu(u) = bu(u) + atom_rates(k);
        bu_support(u) = bu_support(u) + atom_rates_support(k);
    end

    E_users = zeros(U,1);
    for u = 1:U
        idxu = mp_user_atom_indices(atoms, u);
        E_users(u) = sum(P(idxu));
    end
    Et = w(:).' * E_users(:);

    ev = struct();
    ev.bu = bu(:);
    ev.bu_support = bu_support(:);
    ev.atom_rates = atom_rates(:);
    ev.atom_rates_support = atom_rates_support(:);
    ev.E = Et;
    ev.E_users = E_users(:);
    ev.powers = P(:);
    ev.orders = orders;
    ev.active_rx = active_rx;
    ev.rx_rate = rx_rate;
end

%==========================================================================
% ORDER GENERATION / FACE HARVEST
%==========================================================================

function tf = mp_has_tied_theta_cluster(theta, tieTol)
    tf = false;
    U = numel(theta);
    for i = 1:U-1
        for j = i+1:U
            if abs(theta(i)-theta(j)) <= tieTol
                tf = true;
                return;
            end
        end
    end
end

function orderSet = mp_generate_order_candidates(theta, atoms, U, opts, prev_best)
    localOrders = cell(U,1);

    for r = 1:U
        dec = mp_decodable_atoms_at_receiver(atoms, r);
        localW = zeros(numel(dec),1);

        for m = 1:numel(dec)
            k = dec(m);
            u = atoms(k).user;
            base = theta(u);

            if ~isempty(prev_best) && isfield(prev_best, 'rx_rate')
                decset = atoms(k).decoded_at;
                rr = prev_best.rx_rate(k, decset).';
                pos = find(decset == r, 1, 'first');

                if ~isempty(rr) && ~isempty(pos)
                    rmin = min(rr);
                    gap = max(rr(pos) - rmin, 0);
                    closeness = 1 / (opts.orderScoreEps + gap);
                    score = -(base * (1 + opts.orderPrevMix * closeness));
                else
                    score = -base;
                end
            else
                score = -base;
            end

            localW(m) = score;
        end

        localOrders{r} = mp_local_orders_from_scores(dec, localW, opts.orderScoreTol, opts.maxTiePermsPerRx);
    end

    counts = zeros(U,1);
    for r = 1:U
        counts(r) = size(localOrders{r}, 1);
    end
    totalComb = prod(counts);

    if totalComb <= opts.maxGlobalOrderCombos
        orderSet = {struct('rx_orders', {cell(U,1)})};
        orderSet{1}.rx_orders = cell(U,1);

        for r = 1:U
            prevSet = orderSet;
            orderSet = {};

            for a = 1:numel(prevSet)
                for b = 1:size(localOrders{r},1)
                    tmp = prevSet{a};
                    tmp.rx_orders{r} = localOrders{r}(b,:);
                    orderSet{end+1,1} = tmp; %#ok<AGROW>
                end
            end
        end
    else
        nKeep = opts.maxGlobalOrderCombos;
        orderSet = cell(nKeep,1);

        for t = 1:nKeep
            tmp = struct('rx_orders', {cell(U,1)});
            for r = 1:U
                nloc = size(localOrders{r},1);
                if t == 1
                    pick = 1;
                else
                    pick = 1 + mod(t-2, nloc);
                end
                tmp.rx_orders{r} = localOrders{r}(pick,:);
            end
            orderSet{t} = tmp;
        end
    end
        % Add a few random exploratory global order combinations.
    extraSet = cell(opts.nRandomOrderSweeps,1);
    for t = 1:opts.nRandomOrderSweeps
        tmp = struct('rx_orders', {cell(U,1)});
        for r = 1:U
            dec = mp_decodable_atoms_at_receiver(atoms, r);
            tmp.rx_orders{r} = dec(randperm(numel(dec)));
        end
        extraSet{t} = tmp;
    end

    orderSet = [orderSet; extraSet];
end

function O = mp_local_orders_from_scores(atomIdx, scores, tieTol, maxPerms)
    [ws, sidx] = sort(scores(:), 'ascend');
    atomSorted = atomIdx(sidx(:)).';

    groups = cell(0,1);
    gstart = 1;
    for i = 2:numel(ws)
        if abs(ws(i) - ws(i-1)) > tieTol
            groups{end+1,1} = atomSorted(gstart:i-1); %#ok<AGROW>
            gstart = i;
        end
    end
    groups{end+1,1} = atomSorted(gstart:end);

    blockPerms = cell(numel(groups),1);
    for g = 1:numel(groups)
        grp = groups{g};
        if numel(grp) == 1
            blockPerms{g} = grp;
        else
            Pg = perms(grp);
            if size(Pg,1) > maxPerms
                Pg = Pg(1:maxPerms, :);
            end
            blockPerms{g} = Pg;
        end
    end

    O = blockPerms{1};
    for g = 2:numel(blockPerms)
        A1 = O;
        A2 = blockPerms{g};
        newO = zeros(size(A1,1) * size(A2,1), size(A1,2) + size(A2,2));
        tt = 0;
        for i = 1:size(A1,1)
            for j = 1:size(A2,1)
                tt = tt + 1;
                newO(tt,:) = [A1(i,:), A2(j,:)];
            end
        end
        O = newO;
    end

    O = unique(O, 'rows');
    if size(O,1) > maxPerms
        O = O(1:maxPerms, :);
    end
end

function facePool = mp_harvest_tied_theta_face(G, bu_min, w, cb, opts, atoms, theta, prev_best)
    U = numel(theta);

    theta_face = theta(:);
    orderSet = mp_generate_order_candidates(theta_face, atoms, U, opts, prev_best);

    K = numel(atoms);
    facePool = mp_init_vertex_pool(K, U);

    for kk = 1:numel(orderSet)
        cand = mp_inner_solve_for_orders(G, bu_min, w, cb, opts, atoms, orderSet{kk}, theta_face, prev_best);
        facePool = mp_pool_add_vertex(facePool, cand, opts.faceUniqueTol);
    end
end

%==========================================================================
% FINAL GLOBAL SELECTION
%==========================================================================


function [sel, FEAS_FLAG] = mp_global_select(pool, target, hitTol, degTol, opts)
    rawPool = pool;

    % Optional lower-envelope prune for efficiency only.
    % Safety rule: if pruning hurts feasibility, fall back to the raw pool.
    try
        prunedPool = mp_pool_lower_envelope(pool, opts);
        if isempty(prunedPool.bu)
            prunedPool = rawPool;
        end
    catch
        prunedPool = rawPool;
    end

    [sel, FEAS_FLAG] = mp_global_select_core(prunedPool, target, hitTol, degTol, opts);

    if FEAS_FLAG == 0 && size(rawPool.bu,1) > size(prunedPool.bu,1)
        [sel_raw, FEAS_FLAG_raw] = mp_global_select_core(rawPool, target, hitTol, degTol, opts);
        if FEAS_FLAG_raw > FEAS_FLAG || (FEAS_FLAG_raw == FEAS_FLAG && sel_raw.E < sel.E)
            sel = sel_raw;
            FEAS_FLAG = FEAS_FLAG_raw;
        end
    end
end

function [sel, FEAS_FLAG] = mp_global_select_core(pool, target, hitTol, degTol, opts)
    target = target(:);
    U = numel(target);
    K = size(pool.powers, 2);
    nV = size(pool.bu, 1);

    FEAS_FLAG = 0;
    sel = struct();
    sel.mode = 'target-miss';
    sel.bu = zeros(U,1);
    sel.E = inf;
    sel.E_users = zeros(U,1);
    sel.powers = zeros(K,1);
    sel.atom_rates = zeros(K,1);
    sel.atom_rates_support = zeros(K,1);
    sel.best_vertex = struct('idx',[],'bu',[],'E',[],'powers',[],'atom_rates',[],'meta',[]);
    sel.ts_exact = struct('used',0,'kind','','frac',[],'idx',[],'bu_vertices',[],'E_vertices',[],'E_final',[]);
    sel.ts_feasible = struct('used',0,'kind','','frac',[],'idx',[],'bu_vertices',[],'E_vertices',[],'E_final',[]);
    sel.lp_lambda = [];
    sel.lp_active_idx = [];
    sel.lp_active_meta = {};
    sel.lp_exitflag = [];
    sel.lp_message = '';

    if nV == 0
        return;
    end

    bestShort = inf;
    bestIdx = 1;
    for i = 1:nV
        vi = pool.bu(i,:).';
        shorti = sum(max(target - vi, 0).^2);
        if shorti < bestShort || (abs(shorti-bestShort)<=1e-14 && pool.E(i) < pool.E(bestIdx))
            bestShort = shorti;
            bestIdx = i;
        end

        if norm(vi-target, inf) <= hitTol
            FEAS_FLAG = 1;
            sel.mode = 'single-vertex';
            sel.bu = vi;
            sel.E = pool.E(i);
            sel.E_users = pool.E_users(i,:).';
            sel.powers = pool.powers(i,:).';
            sel.atom_rates = pool.atom_rates(i,:).';
            sel.atom_rates_support = pool.atom_rates_support(i,:).';
            sel.best_vertex = mp_pack_best_vertex(pool, i);
            sel.lp_lambda = zeros(nV,1);
            sel.lp_lambda(i) = 1;
            sel.lp_active_idx = i;
            return;
        end
    end

    B = pool.bu;
    Evec = pool.E(:);

    f = Evec;
    Aineq = -B.';
    bineq = -target;
    Aeq = ones(1, nV);
    beq = 1;
    lb = zeros(nV,1);
    ub = ones(nV,1);

    lam = [];
    exitflag = -999;
    lp_message = '';

    try
        if license('test','optimization_toolbox')
            lpOpts = optimoptions('linprog', ...
                'Display', 'none', ...
                'OptimalityTolerance', max(opts.lpTol,1e-9), ...
                'ConstraintTolerance', max(opts.lpTol,1e-9));
            [lam, ~, exitflag, output] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, lpOpts);
            if isfield(output, 'message')
                lp_message = output.message;
            end
        else
            lam = [];
            exitflag = -999;
            lp_message = 'Optimization Toolbox / linprog not available';
        end
    catch ME
        lam = [];
        exitflag = -999;
        lp_message = ME.message;
    end

    sel.lp_exitflag = exitflag;
    sel.lp_message  = lp_message;

    if ~isempty(lam) && exitflag > 0
        lam = mp_compress_lambda_support(lam, B, Evec, target, opts);
        bu_mix = B.' * lam;
        if all(bu_mix >= target - 10*hitTol)
            FEAS_FLAG = 2;
            sel.mode = 'time-share-lp';
            sel.bu = bu_mix;
            sel.E = Evec.' * lam;
            sel.E_users = (pool.E_users.' * lam);
            sel.powers = (pool.powers.' * lam);
            sel.atom_rates = (pool.atom_rates.' * lam);
            sel.atom_rates_support = (pool.atom_rates_support.' * lam);
            sel.best_vertex = mp_pack_best_vertex(pool, bestIdx);
            sel.lp_lambda = lam;
            sel.lp_active_idx = find(lam > max(opts.lpTol,1e-8));
            sel.lp_active_meta = pool.meta(sel.lp_active_idx);

            sel.ts_feasible.used = 1;
            sel.ts_feasible.kind = 'lp-feasible';
            sel.ts_feasible.frac = lam(:);
            sel.ts_feasible.idx = sel.lp_active_idx(:);
            sel.ts_feasible.bu_vertices = pool.bu(sel.lp_active_idx, :);
            sel.ts_feasible.E_vertices = pool.E(sel.lp_active_idx);
            sel.ts_feasible.E_final = sel.E;
            return;
        end
    end

    sel.mode = 'target-miss';
    sel.bu = pool.bu(bestIdx,:).';
    sel.E = pool.E(bestIdx);
    sel.E_users = pool.E_users(bestIdx,:).';
    sel.powers = pool.powers(bestIdx,:).';
    sel.atom_rates = pool.atom_rates(bestIdx,:).';
    sel.atom_rates_support = pool.atom_rates_support(bestIdx,:).';
    sel.best_vertex = mp_pack_best_vertex(pool, bestIdx);
end

%==========================================================================
% MISC HELPERS%==========================================================================
% MISC HELPERS
%==========================================================================

function A = mp_symmetrize_pd(A, theta_max)
    A = 0.5*(A+A.');
    [V,D] = eig(A);
    d = diag(D);
    d = max(d, 1e-14*max(1,max(theta_max.^2)));
    A = V*diag(d)*V.';
    A = 0.5*(A+A.');
end

function out = mp_pack_best_vertex(pool, idx)
    out = struct();
    out.idx = idx;
    out.bu = pool.bu(idx,:).';
    out.E = pool.E(idx);
    out.powers = pool.powers(idx,:).';
    out.atom_rates = pool.atom_rates(idx,:).';
    out.meta = pool.meta{idx};
end

function ts = mp_pack_ts(pool, i, j, alpha, kind)
    ts = struct();
    ts.used = 1;
    ts.kind = kind;
    ts.frac = alpha(:);
    ts.idx = [i; j];
    ts.bu_vertices = [pool.bu(i,:); pool.bu(j,:)];
    ts.E_vertices = [pool.E(i); pool.E(j)];
    ts.E_final = alpha(1)*pool.E(i) + alpha(2)*pool.E(j);
end

function [tf, alpha] = mp_segment_exact_target(v1, v2, t, tol)
    tf = false;
    alpha = [];

    D = v1 - v2;
    idx = find(abs(D) > tol);
    if isempty(idx)
        if norm(v1 - t, inf) <= tol
            tf = true;
            alpha = [1; 0];
        end
        return;
    end

    vals = zeros(numel(idx),1);
    for k = 1:numel(idx)
        ii = idx(k);
        vals(k) = (t(ii) - v2(ii)) / D(ii);
    end

    a1 = vals(1);
    if max(abs(vals - a1)) > 1e-5
        return;
    end

    a2 = 1 - a1;
    mix = a1*v1 + a2*v2;

    if a1 >= -tol && a2 >= -tol && norm(mix - t, inf) <= 10*tol
        tf = true;
        alpha = [a1; a2];
    end
end

function [tf, lamLo, lamHi] = mp_segment_dominates_target(v1, v2, t, tol)
    U = numel(t);
    lamLo = 0;
    lamHi = 1;

    for k = 1:U
        d = v1(k) - v2(k);
        if abs(d) <= tol
            if v2(k) < t(k) - tol
                tf = false;
                lamLo = NaN; lamHi = NaN;
                return;
            end
        elseif d > 0
            lamLo = max(lamLo, (t(k)-v2(k))/d);
        else
            lamHi = min(lamHi, (t(k)-v2(k))/d);
        end
    end

    lamLo = max(lamLo,0);
    lamHi = min(lamHi,1);
    tf = (lamLo <= lamHi + tol);
end

function lam = mp_closest_lambda_to_target(v1, v2, t, lamLo, lamHi)
    d = v1-v2;
    if norm(d) <= 0
        lam = lamLo;
        return;
    end
    lam0 = (d.'*(t-v2)) / (d.'*d);
    lam = min(max(lam0, lamLo), lamHi);
end


function lam = mp_compress_lambda_support(lam, B, Evec, target, opts)
    lam = lam(:);
    if ~isfield(opts,'compressActiveSet') || ~opts.compressActiveSet
        return;
    end
    if ~license('test','optimization_toolbox')
        return;
    end

    tol = max(opts.lpTol, 1e-9);
    active = find(lam > max(tol,1e-8));
    U = size(B,2);

    % On a proper supporting face in R^U, at most U active points are needed.
    % We greedily test whether one active point can be removed while preserving
    % feasibility and objective value.
    while numel(active) > U
        dropped = false;
        for kk = 1:numel(active)
            idx = active;
            idx(kk) = [];
            Bk = B(idx,:);
            Ek = Evec(idx);

            f = Ek(:);
            Aineq = -Bk.';
            bineq = -target(:);
            Aeq = ones(1, numel(idx));
            beq = 1;
            lb = zeros(numel(idx),1);
            ub = ones(numel(idx),1);

            try
                lpOpts = optimoptions('linprog', ...
                    'Display','none', ...
                    'OptimalityTolerance', tol, ...
                    'ConstraintTolerance', tol);
                [lamk, valk, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, lpOpts);
                if exitflag > 0
                    buk = Bk.' * lamk;
                    if all(buk >= target(:) - 10*tol) && valk <= Evec(:).' * lam + 10*tol
                        lamNew = zeros(size(lam));
                        lamNew(idx) = lamk;
                        lam = lamNew;
                        active = find(lam > max(tol,1e-8));
                        dropped = true;
                        break;
                    end
                end
            catch
            end
        end
        if ~dropped
            break;
        end
    end
end

function y = mp_tern(cond, a, b)
    if cond
        y = a;
    else
        y = b;
    end
end

function sched = mp_build_schedule_from_selection(sel, ~, ~, ~)

    if isfield(sel,'lp_active_meta') && ~isempty(sel.lp_active_meta)
        nS = numel(sel.lp_active_meta);
        idx = sel.lp_active_idx(:);
        lam = sel.lp_lambda(idx);
    else
        nS = 1;
        lam = 1;
    end

    sched = repmat(struct( ...
        'lambda', [], ...
        'pool_index', [], ...
        'powers', [], ...
        'atom_energies', [], ...
        'atom_rates', [], ...
        'atom_rates_support', [], ...
        'bu', [], ...
        'E', [], ...
        'E_users', [], ...
        'orders', [], ...
        'rx_orders', [], ...
        'rx_rate', [], ...
        'active_rx', []), nS, 1);
 
    for s = 1:nS
        if isfield(sel,'lp_active_meta') && ~isempty(sel.lp_active_meta)
            meta = sel.lp_active_meta{s};
            sched(s).lambda = lam(s);
            sched(s).pool_index = idx(s);
        else
            meta = sel.best_vertex.meta;
            sched(s).lambda = 1;
            sched(s).pool_index = sel.best_vertex.idx;
        end

        sched(s).powers = meta.powers(:);
        sched(s).atom_energies = meta.powers(:);
        sched(s).atom_rates = meta.atom_rates(:);
        sched(s).atom_rates_support = meta.atom_rates_support(:);

        sched(s).bu = meta.bu(:);
        sched(s).E = meta.E;
        sched(s).E_users = meta.E_users(:);

        sched(s).orders = meta.orders;
        sched(s).rx_orders = meta.orders.rx_orders;
        sched(s).rx_rate = meta.rx_rate;
        sched(s).active_rx = meta.active_rx;
    end
end